Tunneling through an Anderson impurity between superconductors 
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We consider an Anderson impurity (A) weakly connected to a superconducting electrode (S) on 
one side and a superconducting or a normal metal electrode (A*') on the other side. A general path 
integral formalism is developed and the response of SAN and SAS junctions to a constant voltage 
bias V is elucidated, using a combination of the Keldysh technique (to handle non- equilibrium 
effects) and a dynamical mean field approximation (to handle repulsive Hubbard interactions). An 
interesting physics is exposed at sub-gap voltages {eV < A for SAN and eV < 2 A for SAS). For an 
SAN junction, Andreev reflection is strongly affected by Coulomb interaction. For superconductors 
with p-wave symmetry the junction conductance exhibits a remarkable peak at eV < A, while for 
superconductors with s-wave symmetric pair potential the peak is shifted towards the gap edge 
eV = A and strongly suppressed if the Hubbard repulsive interaction increases. Electron transport 
in SAS junctions is determined by an interplay between multiple Andreev reflection (MAR) and 
Coulomb effects. For s-wave superconductors the usual peaks in the conductance that originate from 
MAR are shifted by interaction to larger values of V. They are also suppressed as the Hubbard 
interaction strength grows. For p-wave superconductors the sub-gap current is much larger and the 
I — V characteristics reveal a new feature, namely, a peak in the current resulting from a mid-gap 
bound state in the junction. 
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I. INTRODUCTION 

The dynamical behavior of Josephson junction 
strongly depends, among other factors, on its trans- 
parency. If the insulating barrier is not too high then the 
concept of nonlinear tunneling becomes relevant. In this 
case the characteristic dynamical conductance dl/dV at 
applied voltages V less than the superconducting gap 
A shows a sub-gap structure. An explanation of this 
behavior was given some time ago based on the 

mechanisms of multiple Andreev reflections (MAR). Re- 
cently, the sub-gap current was calculated for the case 
of electron tunneling through a junction with resonant 
impurity Rapid progress in the technology of su- 
perconducting junctions makes it possible to fabricate 
junctions composed of quantum dots weakly coupled to 
superconducting or normal electrodes. The basic physics 
of such a device can be elucidated once it is modeled as 
an Anderson impurity center. In this case the Coulomb 
interaction is expected to strongly affect the tunneling 
current in general and the sub-gap current in particu- 
lar. Since the sub-gap current is originated from multi- 
ple Andreev reflections, its physics has a close similar- 
ity to that of the Josephson current. In this context, 
it is established |^ that the tunneling through a quan- 
tum dot is suppressed if the effective Kondo temperature 
Tk =VUT exp [— 7r|eo|/2r] is small as compared to the 
superconducting gap A. Hereafter, U is the Hubbard 



repulsion strength, eq is the orbital energy of the dot 
electron and T is the width of this energy state. Strong 
interaction-induced suppression of the current through 
superconducting quantum dots was also observed exper- 
imentally j5|. 

Quite recently detailed measurements of the I — V 
curves in atomic-size metallic contacts were performed 
Q . An explanation of the observed I — V curves were 
given [0 in terms of the atomic valence orbitals which 
represent different conducting channels. The Coulomb 
interaction was considered there to be screened as in bulk 
metals. However, for quantum dots and break junctions 
the screening is virtually ineffective and an unscreened 
Hubbard-type repulsive interaction emerges. In this case 
the Kondo temperature Tff becomes a relevant parame- 
ter, separating levels with > A (which are responsible 
for high, nearly resonant conductance) from levels with 
Tk < A, in which the conductance is strongly influenced 
by interaction. 

One of the main goals of the present paper is to de- 
velop a detailed theoretical analysis of an interplay be- 
tween the phenomena of multiple Andreev reflections and 
Coulomb interaction in superconducting quantum dots. 
Even though both MAR and Coulomb effects have been 
intensively studied in the literature over past decades, 
an interplay between them - to the best of our knowl- 
edge - was not elucidated until now. In this paper we 
will demonstrate that these two phenomena, being com- 
bined in superconducting quantum dots, lead to novel 
physical effects and - depending on parameters - may 
dramatically inffuence the sub-gap conductance pattern 
of the system. In short. Coulomb suppression of MAR 
turns out to be much more pronounced than, say, that of 
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single electron tunneling. This is because during MAR 
cycles at relatively low voltages the charge is transferred 
by large quanta, much larger than the electron charge e. 

Another important issue in the study of SAS and SAN 
junctions is the parity of the order parameter of the su- 
perconducting electrodes. For example, the order param- 
eter of the recently discovered ^ superconducting mate- 
rial Sr2RuOi is believed to have a p- wave symmetry 
If a superconductor of this type is properly oriented with 
respect to the tunneling direction the principal contribu- 
tion to the Josephson current comes from a bound state 
pO| , pl| formed at the contact point. This bound state 
arises since the pair potential has an opposite sign for in- 
jected and reflected quasiparticles and is expected to play 
an important role in the formation of sub-gap currents. 

In the present work we expose the physics of SAS and 
SAN junctions subject to a finite potential bias. In par- 
ticular, we calculate the tunneling current and the dy- 
namical conductance for junctions consisting of s- and 
p-wave superconductors. The main steps required for 
treating the pertinent many-body problem can be sum- 
marized as follows: 1) Taking the Fermi energy of the 
unbiased lead as an energy reference, the site energy eo 
of the Anderson impurity is chosen such that eo < while 
U + €() > 0. These inequalities assure that assuming the 
quantum dot to be at most singly occupied should be an 
excellent approximation. 2) To handle the strong inter- 
action appearing in the Hubbard term, a mean field ap- 
proximation 1 12 13] is adopted. 3) The formalism should 



take into account the nonequilibrium nature of the phys- 
ical system. For this purpose, the standard approach is 
to start from the expression for the kernel of the evolu- 
tion operator or the generating functional, which is the 
analog of the partition function in the equilibrium case, 
evaluated, however, on a Keldysh contour ||l^ (see review 
article in Ref. |15|). At the end of this procedure one is 
able to calculate the SAN Andreev conductance analyt- 
ically, and to get expressions for the non-linear response 
of SAS junctions which are amenable for numerical eval- 
uation. 

The technical procedure by which we manage to ad- 
vance the calculations is detailed below in section 2, 
where we derive an effective action for SAS and SAN 
junctions. In section 3 we discuss the dynamical mean 
field approximation adopted in the present work in order 
to treat interaction effects. Concrete results pertaining 
to sub-gap current in SAS junctions and differential con- 
ductance in SAN junctions are presented and discussed 
in section 4. The paper is then concluded and summa- 
rized in section 5. Some technical details of the calcula- 
tion are given in Appendix. 



A. The Model 

Consider a system consisting of two superconducting 
wide strips on the left [x < 0, — oo < y < oo) and on 
the right {x > 0, — oo < y < oo) weakly connected by a 
quantum dot through which an electron tunneling takes 
place. This system can be described by the Hamiltonian 



H = H, 



Hp + H 



dot 



(1) 



The Hamiltonians of the left and right superconducting 
electrodes have the standard BCS form 



H, = I dr[vI/t^(r)e(V)*,,(r) 



(2) 



Here ^E*]^ i'^ja) are the electron creation (annihilation) 
operators, A is the BCS coupling constant, ^(V) = 
— V^/2m — /i and j — L,R. Here and below we set 
the Planck's constant h = 1. Whenever appropriate, the 
spin, space and time dependence of all the field operators 
will not be explicitly displayed. 

The quantum dot is treated as an Anderson impurity 
center located at x = y ~ 0. It is described by the 
Hamiltonian 



H 



dot 



(3) 



where C| and Co- are the electron operators in the dot. 
The impurity site energy eo (counted from the Fermi en- 
ergy /x) is assumed to be far below the Fermi level eo < 0. 
The presence of a strong Coulomb repulsion J7 > — eo be- 
tween electrons in the same orbital guarantees that the 
dot is at most singly occupied. 

Electron tunneling through the dot is accounted for by 
means of the term, 



(4) 



J=L,R 



where Tl(r) are the effective transfer amplitudes between 
the left (right) electrode and the dot. 

In what follows we will always assume that, if a bias 
voltage V is applied to the system from, say, right to left, 
the entire voltage drop occurs across the dot. Hence, the 
quasiparticle distribution functions in the leads are the 
Fermi ones, with the chemical potentials of the electrodes 
shifted with respect to each other by eV . 



B. Evolution Operator 



II. GENERAL ANALYSIS 



Complete information about the quantum dynamics of 
the system is contained within the evolution operator de- 
fined on the Keldysh contour [H K (which consists of 
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forward and backward oriented time branches). The ker- 
nel J of this evolution operator can be expressed in terms 
of a path integral, 



'dot 



dt 



dt 



U 



{ccf 



(11) 



(5) 



So^ [dt f dr^,{r,t)Gj'^j{r,t) 



over the fermion fields corresponding to the operators 
and C (here the field ^ corresponds to 



j=L.R 



LV ii' -RT' Rl 



) and similarly for other fields), S 



(12) 



Jj^ Ldt is the action and L is the Lagrangian pertaining Here we defined e ^ eg + U/2. In order to obtain the 

expression for the operator we employ the standard 



to the Hamiltonian (|l|). The external fields (e.g. electro- 
magnetic fields) can be treated as the source terms for 
the action, though, the fluctuating parts of these fields 
should be integrated as well. 

Usually it is convenient to perform an operator rota- 
tion C c and — > -0 in Keldysh space: 



Ca,Q- 



QC; ij = ^a,Q'\ V = Q*- (6) 



Here cr^ is one of the Pauli matrices (Jy, operating 
in the Keldysh space and 



(7) 



is the Keldysh matrix. The Grassman variables c, c, if), 
■0 are now defined solely on the forward time branch. 

The transformation of the Green functions follows di- 
rectly from (j|). One starts from the 2x2 matrix G of 
the Green functions defined in terms of the initial electron 
operators. The elements of the matrix G are the Green 

functions Gij with i,j = +, — according to whether 
the time belongs to the upper or the lower branch of 
the Keldysh contour K. Of these four Green functions 
only three are independent. Under the operator rota- 
tion (^) the Green-Keldysh matrix G is transformed as 
G = Q-^GQ, where 



G 



[ 6 G^ 



(8) 



and 



G« = -ie{t-t') {i}{r,t)i;\r' ,t')+il;\r\t')i}{r,t)) , 
G^ = i9{t' - t) (0(r, t)V/(r', t') + V/(r', t')V(r, t)) , 
= (^(r, i)V'^(r-', t') + V'^(r', t')V'(r, t)) , (9) 

are respectively retarded, advanced and Keldysh Green 
functions. Each of these matrices is in turn 2x2 matrix 
in the Nambu space. 

The path integral ^ is now expressed in terms of the 
new Grassman variables 

J = J V'^V'ipVcDcexp{iSdot+iSo[^,^/j]), (10) 
where 



Hubbard-Stratonovich transformation of the ^'^-term in 
(^) and introduce additional path integrals over the com- 
plex scalar order parameter field A{r,t) defined on the 
Keldysh contour, see e.g. Ref. Since here we are not 
interested in the fluctuation effects for the order param- 
eter field, we will evaluate the path integral over A by 
means of the saddle point approximation, which amounts 
to setting A{r,t) equal to the equilibrium superconduct- 
ing order parameter values A^n of the left and the right 
electrodes. If needed, fluctuations of the order parame- 
ter field (both the amplitude and the phase) can easily be 
included into our consideration along the same lines as 
it was done in Ref. Disregarding such fluctuations 

here, we flnd 



Gz^fl(0 = »|-T.e(v) 



L.B. 



t^^Ib. (13) 



where we define t± — [t^ ± iTy)/2. Here and below 
Tx,Ty,Tz is the set of Pauli matrices operating in the 
Nambu space (for the sake of clarity we chose a different 
notation from that used for Pauli matrices operating in 
the Keldysh space). 



C. Effective Action 

Let us now proceed with the derivation of the effective 
action for our model. We first notice that the ^-fields de- 
pendent part 5*0 of the total action is quadratic in these 
fields. Hence, the integrals over ip and ip in ( |l0|) can be 
evaluated exactly, resulting in an action S'cnv(c, c) defined 



exp(iS'o„v[c, c]) = / 'D^p'Dil;exp{iSo[ilj,ip]). 



(14) 



Its physical content can be understood as follows; One 
can say that electrons in two superconducting bulks serve 
as an effective environment for the quantum dot. Inte- 
grating out these electron variables in the spirit of the 
Feynman- Vernon influence functional approach [ p^ one 
arrives at the "environment" contribution to the action 
Scnv expressed only in terms of the Anderson impurity 
variables c and c. 
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Due to the fact that coupHng to the leads is concen- 
trated at one point {x,y) = (0,0) we can integrate out 
the fields inside the superconductors (hereafter referred 
as bulk fields) and obtain an effective action in terms 
of fermion operators with arguments solely on the sur- 
face. In order to achieve this central goal let us first note 
that translation invariance along y permits the Fourier- 
transform in eq. ( p^ ) in this direction. The problem 
then reduces to a one dimensional one with fermion fields 
ipkix) where k is the momentum along y. Gaussian inte- 
gration over the bulk fields can be done with the help of 
the saddle point method. 

Let us consider, say, the left superconductor and omit 
the subscript j = L for the moment. The pertinent equa- 
tion for the optimal field reads, 



r(e) = (.r(e)-r(e))tanh(e/2r). 



(21) 



(15) 



where £,x — — (l/2m)(9^/5a;^) — /i/t and 11^=11— k^ /2m. 

Let us decompose ipk{x) — i^k^x) + -0^(0) in such a 
way that on the surface one has '0^(0) = 0- The bulk 
field '4'\{x) satisfies the inhomogeneous equation: 

G-i(ex)0fc(x) = -AifcT,0fc(O). (16) 

In the right hand side of this equation we employed the 
standard quasiclassical (Andreev) approximation which 
makes use of the fact that the superconducting gap as 
well as other typical energies of the problem are all much 
smaller than the Fermi energy. 

In order to solve eq. ( p^ we find the Green function 
Gk{x,t; x' ,t') (which satisfies the same equation albeit 
with a S -function on the right hand side) and require 
G to vanish at x = 0. The solution of eq. (|l^) is then 
exploited to express ip'lix) in terms of the surface fields 



■0^(0). Combining the result with eq. ( p2[ ) we arrive at 
the intermediate effective action S for a superconducting 
electrode which depends on the ■0-fields at the surface, 

S^i j dt j dt'Y, y ^fe(0, t)T,9{t, t')MO, t'), (17) 



Here ip{t) — (po + 2e J* V{ti)dti is the time-dependent 
phase of the superconducting order parameter and V{t) 
is the electric potential of the superconducting electrode. 

An identical procedure applies for the right elec- 
trode. Each superconductor is thus described by a zero- 
dimensional action, respectively Sl and 5*^, coupled by 
an on-site hopping term with the Anderson impurity. It 
is now possible to integrate out these surface fields. The 
integral 



(22) 



can easily be evaluated, so that the contribution of the 
superconductors to the total effective action of our model 
is manifested in Scnv defined as 

S'cnv ^2iJdtJ dt'J2cit)T,[^gL{t, t') 



J = J P0(O)I?V(O) eM^SL,R 
+i / dt{TLMMO)T.c + c.c.)) 



-Hgj,{t,t')]c{t'). 



(23) 



Note that in deriving (^3|) we made use of the normaliza- 
tion condition ll^ 9l r — 1- 

Eq. ( p^ ) is valid for an arbitrary pairing symmetry. 
In the case of unconventional superconductors the Green 
functions gL,R. depend explicitly on the direction of the 
Fermi velocity. For uniform s-wave superconductors such 
dependence is absent and eq. ( |2^ ) can be simplified fur- 
ther. Defining the tunneling rates between the left (right) 
superconductor and the dot as 



T2 

r 1MB. 

k 



(24) 



where Vx = •\/2/ifc/m is the quasiparticle velocity in the 
x-direction. For a uniform superconducting half-space 
(here the left one), the Grecn-Kcldysh matrix 

9it^t')T, = -i^— J dx'Gk{x,t;x',t')\x=o (18) 

(which has the structure (^)) is expressed in terms of the 
Eilcnberger functions (EtI as follows: 



we obtain 



de ifit')^:, 

e 2 , (19) 



where 



gR/Af^^^^ _ (e±iO)r^-)-i|A|Ty 



v/(e±*0)2-|A|2' 



(20) 



5cnv = :^ / dt / dt'c{t)T,[TL9L{t.t') + TRgR{t,t')\c{t'). (25) 



The definition ( p4| ) requires a comment. Here we consider 
the situation with one conducting channel in the dot. In 
this case the transfer amplitudes Tl^r should effectively 
differ from zero only for \vx \ ~ vp. One can easily gener- 
alize the action (^5|) to the situation with several or even 
many conducting channels. In this case the summation 
over momentum (essentially equivalent to the summa- 
tion over conducting modes) should be done in ( p5| ) and 
some other dependence of ^ on Vx should apply. For 
instance, for tunnel junctions in the many channel limit 



one can demonstrate [|l8| that Tl 



It is also quite 



4 



clear that the transfer amplitudes T^^r cannot be con- 
sidered as constants independent of the Fermi velocity 
direction, as it is sometimes assumed in the literature. 



In that case the sum (24) would simply diverge at small 
Vx in a clear contradiction with the fact that quasipar- 
ticles with Vx ^ ^ should not contribute to the current 
at all. This "paradox" is resolved in a trivial way: the 
amplitudes T^ r do depend on Vx and, moreover, they 
should vanish at w^; — + 0. For further discussion of this 
point we refer the reader to Ref . . 

Combining eqs. ( p^ ) and (|l^) we arrive at the expres- 
sion for the kernel of the evolution operator J solely in 
terms of the fields c and c: 



J = / I?cI?cexp(zS'cff), 5'off[c, c] = S'dot + Sc 



(26) 



Here S'cff[c, c] (defined by eqs. ( [Tl| ) and (|25|)) represents 
the effective action for a quantum dot between two su- 
perconductors. 



where the path integral J is defined in (22). The func- 
tional derivative of (p9| ) with respect to the ry-fields just 
yields the function G^: 



Gni 



. 5^Z 
SfjSri 



\ri=ri=0- 



(30) 



Evaluating the path integral (|29|) and making use of (|3C 
we arrive at the following identity 



iG^ = 



-QLTz 



(31) 



Combining (28) and ( pl| ) with the condition g\ = 1 
we observe that the contribution of the first term in the 
right-hand side of ( pT| ) to the current vanishes identically, 
and only the second term oc (cc) turns out to be impor- 
tant. Making use of the definition ( [2^ ) and symmetrizing 
the final result with respect to R and L we arrive at the 
following expression for the current 



D. Transport current 



-TriiTLgL ~^r9r){cc)\k +h.c.]. 



(32) 



In order to complete our general analysis let us ex- 
press the current through the dot in terms of the correla- 
tion function for the variables c and c. This goal can be 
achieved by various means. For instance, one can treat 
the superconducting phase difference across the dot as a 
source field in the effective action and obtain the expres- 
sion for the current just by varying the corresponding 
generating functional with respect to this phase differ- 
ence. Another possible procedure is to directly employ 
the general expression for the current in terms of the 
Green-Keldysh functions of one (e.g. the left) supercon- 
ductor, with arguments at the impurity site: 



1 = 



Am 



dy{dx ~ dx')'Tr[G'' {xy,x'y'-t)]x 



(27) 



where the trace is taken in the Nambu space. 

As before, it is convenient to split the bulk and the 
surface variables. After a simple algebra we transform 
eq. (p7|) into the following result: 



I = -ij^VxTT[gLG^ - h.c.]|if , 



(28) 



where G^ = —i < ipk{0)'ipk{0) > is the Green-Keldysh 
function for the surface ?/;-fields. Here and below the 
integration over the internal time variables in the prod- 
uct of matrices is implied and {■■■)\k means the Keldysh 
component of this product. 

Finally, let us express the function G^ in terms of the 
correlator for the fields c and c. Let us consider the gen- 
erating functional for the surface fields 



Z[ri,r]] = J[TlTzC + ti^Tlt^c + ri], 



(29) 



This expression completes our derivation. We have 
demonstrated that in order to calculate the current 
through an interacting quantum dot between two super- 
conducting electrodes it is sufficient to evaluate the cor- 
relator (cc) in the model defined by the effective action 
S'cff = S'dot + (prl), (p5|). Our approach enables one 
to investigate both equilibrium and nonequilibrium elec- 
tron transport in superconducting quantum dots. In the 
noninteracting limit J7 — > the problem reduces to a 
Gaussian one. In this case it can easily be solved and, 
as we will demonstrate below, the well known results de- 
scribing normal and superconducting contacts without 
interaction can be re derived in a straightforward man- 
ner. In the interacting case U ^ the solution of the 
problem involve approximations. One of them, the dy- 
namical mean field approximation, is described in the 
next section. 



III. MEAN FIELD APPROXIMATION 

In order to proceed further let us decouple the inter- 
acting term in (O) by means of a Hubbard-Stratonovich 
transformation |n2 13 1 introducing additional scalar fields 
7-1- . The kernel J now reads, 



J 



I?cl?cl?7+I?7_ exp 



dtc i 



dt 



(33) 



dt ( c7+cr^c + c^^c - 



(34) 
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These equations are still exact. Now let us as- 
sume that the effective Kondo temperature Tk 
=y/UT exp [— 7r|eo|/2r] is smaller than the superconduct- 
ing gap A. In this case interactions can be accounted for 
within the mean field (MF) approximation. The fields 
7± in (jS^) can be determined from the saddle point con- 
ditions 



5J/5-f± = 0. 



(35) 



In general these two equations contain an explicit depen- 
dence on the time variable. Let us average these equa- 
tions over time and consider 7-1- as time independent pa- 
rameters. This approximation is equivalent to retaining 
only the first moment of 7±. The self-consistency equa- 
tions (^) now read 



7+ = ^ / at <cc>, 



(36) 
(37) 



As it turns out from our numerical analysis (to be de- 
scribed below), the parameter 7+ has a negligible ef- 
fect on the sub-gap current. It just slightly renormal- 
izes the coupling constants Tl^r of our model. On the 
other hand, the second parameter, 7_, has a large im- 
pact on the I — V characteristics. Therefore in what 
follows we will set 7+ = and take into account only the 
second self-consistency equation { ^t\ ) for 7_. Under this 
approximation the effective action of our model acquires 
the following form 

de 



Soffb] ^ J de'cM{e,e')c, 

M(e, e') - ^ (e - e') [e + 7- - + ir, iTji/2)gnie)] 
+tT,{TL/2)gL{e,e'). (38) 

Here and below we deliberately choose the electrostatic 
potential of the right electrode to be equal to zero. In 
this case the Keldysh matrix gn is diagonal in the energy 
space. Performing the functional integration over Grass- 
man variables c and c we can cast the self-consistency 
equation ( ^7|) for 7_ in terms of the matrix 



(M 



where M^, are three independent elements of 

the Keldysh matrix M (^8|). Recall that each of these 
elements is a 2 x 2 matrix in the Nambu space and an 
infinite matrix in the energy space. Eq. (|37| ) for 7_ can 
now be rewritten as 



7_ = i^TriM'')- M {Mn~\ 



(40) 



with the trace being taken both in the energy and spin 
spaces. 



Finally, employing the MF approximation for the Hub- 
bard interaction as was implied in the calculation of 7_ , 
we get the current as a difference of symmetric forms. 



(41) 



Consider now the case of a constant in time voltage 
bias V and recall that the entire voltage drop occurs 
across the quantum dot. Setting the phase of the right 
electrode equal to zero, for the phase of the left super- 
conductor we obtain cpit) — 2eVt + ipQ. Let us express 
gL in terms of the matrix elements in the energy space 



{e\gLW)= J2 S{e~e' + 2seV)gL{e,e + 2seV), 

s=0,±l 

gi(e, e + 2seV) = {gi\E - eV)P+ + gf{e + eV)P^)So,s 
+e''^"giHe ~ eV)T+6s^^, + e-'^'gl^e + eV)T_6ss, (42) 

where the superscripts denote the matrix elements in the 
Nambu space and P± = (1 ± r2)/2. 

In what follows we shall abbreviate gL(e -I- 2rneV, e -f 
2neV) = (m|gi(e)|n) where the right hand side is ob- 
tained from eq. ( |42[ ) after replacing e — > e 4- 2meV, 

1. Then 



we have 



-'n,m-\- 



gLie, e') = ^5(e - e' + 2neV)iO\gL{e)\n), (43) 

n 

The matrix M (|38| ) may also be represented in a sim- 
ilar form, that is, 

M(e,e') ==^d{e - e' + n2eV){0\M{e)\n), (44) 

n 

where 

(m|Af(e)|n) = 6„i.n[e + m2eV + 7^ 



'Tzi - ^^TzgB.{e + Tn2eV)] 



\m\gL{e)\n). 



(45) 



The integration over energy variables in the self- 
consistent equation for 7_ and in the expression for 
the time averaged current is conveniently performed by 
dividing the whole energy domain into slices of width 
2eV and performing energy integration on an interval 
[0 < -E < 2eV] . Thus we can use the discrete representa- 
tion (Ba) and write 



7- 



4 r ^Y.^<^\^^i'')-'M^{M^r'\n). 



8 Jo 27r ^ 



-{L^ R)])+h.c.]\n). 



(46) 
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Let us also note that in the case of SAN junctions the 
expressions for the current and for 7_ can be simphfied 
further. E.g. it is easy to observe that in this case eq. 
( p] ) takes the form 



X= U + 7- + 



2 



(51) 



In the hmit eT^ < A and T 
G = I/V we obtain 



de. 



-Tr[((M«(e))-\/(6,y)(Af4(6))-\gfC 



2tt 



G 



where the matrix / has the standard form 

, /tanh(^) 

J^'^^>-[ tanh(^) 



(47) 



h /rHr? 



for the conductance 



(52) 



+ 7-- 



(48) 



Eq. (^^ can be straightforwardly evaluated since the 
(Fourier transformed) matrices {M^'^)~^ depend now 
only on one energy e {g^ in (|38| ) is proportional to 6{e—e') 
in this case) and, hence, can easily be inverted analyti- 
cally. Similar simplifications can also be performed in 
the self-consistency eq. (^0[). 



IV. RESULTS AND DISCUSSION 
A. The SAN junction 

1. s-wave superconductors 

We start from calculating the differential conductance 
of an SAN contact assuming the s-wave pairing symme- 
try in a superconducting electrode. As it was already 
pointed out above, eqs. (^), (48) allow one to proceed 
analytically. From these equations one obtains the ex- 
pression for current which consists of two parts. The first 
part originates from the integration over sub-gap energies 
e < A and yields the dominating contribution to the cur- 
rent at low temperatures. The other part comes from 
integration over energies e > A. At low voltages and 
temperatures (lower than the gap A) this second part 
gives a negligible contribution to the current. Consider- 
ing below the sub-gap contribution only, we find 



eFLFij 



tanh 



eV\ , fe-eV 

- tanh I 



2T 



(49) 



where at sub-gap voltages and energies one has 
AM|A|-k|) 



— 



El 

4 



A2 



4 



(50) 



and 



In order to recover the expression for G in the noninter- 
acting limit in eq. ( ^^ one should simply put 7_ = 0. 
In a symmetric case Tl — and for e eq. ( [s^ ) 
reduces to the well known result ||l^ 



Gms — 2Gnn — —r- 



4e^ 
~h 



(53) 



In the presence of Coulomb interaction the parame- 
ter 7_ in ( po[^ ) should be determined from the self- 
consistency equation (^0[). This equation has a solution 
provided the interaction U is not very large, i.e. outside 
the Kondo regime. In the limit of large U the solution of 
eq. ( ^0| ) is absent, which indicates the failure of the MF 
approximation in the Kondo regime. 

Here eq. (^0|) was solved numerically for a given set 
of the system parameters. In our numerical analysis we 
chose Tl = F^j = F = 0.35A and considered the most 
interesting sub-gap voltage bias regime eV < 2A in the 
low temperature limit T —^ 0. The values of the Hub- 
bard repulsion parameter U were fixed to he U — 2.45A 
and U = 2.72A. For convenience we scaled our equations 
expressing the parameters e, U, h and T in units of 
A. After that the current and the conductance are nor- 
malized respectively in units of Ae/h and /2h. 

For the above parameters the value 7_ was found to 
be 7_ « 1.2 A with variations (depending on the voltage) 
within 10 15 per cents. Thus within a reasonable ac- 
curacy in eqs. (p0| - |52|) one can consider 7_ as constant. 
On the other hand, even though the variations of 7_ 
with voltage are not large in magnitude, sometimes they 
occur over a small voltage interval. Therefore such varia- 
tions may have a considerable impact on the differential 
conductance cr = dl/dV and should also be taken into 
account. This was done within our numerical analysis. 
The corresponding results are presented in figure 0. 

We observe that for given parameters the conductance 
virtually vanishes in the substantial part of the sub-gap 
region. Note, however, that at voltages close to but still 
sma/Zer than A/e the differential conductance a increases 
sharply. This feature can easily be understood as a result 
of interplay between Coulomb blockade and two electron 
tunneling effects. It is well known that the sub-gap 
conductance in SN junctions is caused by the mecha- 
nism of Andreev reflection during which the charge 2e is 
transferred between the electrodes. Without interaction 
eq. (|5^) holds down to V —i- while in the presence of 
interaction and at T = two electron tunneling process 
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□ 



□ 



FIG. 1. Differential conductance of an SAN junction with 
s-wave symmetry superconductor. The figure displays the 
dependence of Andreev conductance on the applied voltage 
(in units of superconducting pair potential) for U = 2.450 
(solid line curve G) and U — 2.713 (dotted line curve D). The 
barrier transparency is F = 0.35 and the dot level energy is 
eo = -1.5. 



A similar study of an interplay between two-electron 
tunneling and Coulomb effects in SIS junctions was car- 
ried out in Ref. pi. 



There exists also a certain analogy between our results 
and those obtained for superconductor-ferromagnet (SF) 
junctions ||2^]. Here the repulsion parameter U plays a 
role similar to that of an exchange term in SF systems: 
in both cases the sub-gap conductance can be tuned by 
changing this parameter in a way that a smaller value of 
U corresponds to a weaker exchange field. In contrast to 
our system, however, changing of the exchange field in 
SF junctions leads to smooth variations of the sub-gap 
conductance [Esl. 



Let us also note that here we do not consider the Kondo 
limit ]26|-p^ , in which case a zero bias conductance peak 
is expected. This peak appears, simply, because in the 
Kondo limit and at = there exists an open channel 
between the normal metal and the superconductor. Here 
the Kondo temperature Tk is assumed to be small and, 
hence, the zero bias peak is absent. 

Let us now briefly consider the limit of large bias volt- 
ages eV ^ A. In this case the current may be repre- 
sented as a sum of two terms I = I\+ 12- The term Ii is 



determined by the expression similar to (49) which now 



includes the contribution from energies above the gap. 
We find 



h = 
4 

tanh 



f°° de 

J-oc 



eV 



2T 



tanh 1^ 



e-eV 



2T 



(54) 



Here B{e) is again given by eq. (|5 
Bi{e) reads 



while the function 



2m\e\~\A\) 



+ (e + 7-)'+Xi] 



i^' + 1-? + Xi f + + l-nVL + TR^^y ' 



(55) 



Here we also defined 
Xi 



(56) 



The other contribution I2 is proportional to the level po- 
sition e. One obtains 

,2^ 



I2 = 



62(e)[tanh(^i^) + tanh(i^) - 2 tanh(^)]. (57) 

The expression for B2{e) can be obtained from eq. ( [55| ) 
if one replaces the term in the square brackets by the 
expression — 2e(e + 7-). 

The above results together with the self-consistency 
equation for 7_ provide a complete description for the 
I —V curve of an SAN junction in the presence of inter- 
actions. In all interesting limits the energy integrals in 
(^), (^) can be carried out and the corresponding ex- 
pressions for the current can be obtained. These general 
expressions, however, turn out to be quite complicated 
and will not be analyzed in details further below. 

Here we just demonstrate that in the noninteracting 
limit 7_ = our results reduce to the results already 
well known in the literature. In the leading order approx- 
imation eqs. (54), ( ^5| ) yield the standard Breit-Wigner 
formula 



2e^ 



TlTr, 



(58) 
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FIG. 3. Same as in figure ^ but for a p-wave symmetry su- 
perconductor. The parameters U, T and eo are identical to 
those of figure nl 



FIG. 4. The subgap tunneling current versus voltage for an 
SAS junction for which the superconductors pair potential 
has an s-wave symmetry. The parameters are U = 2.4, 2.7, 
e = -1.5 and F = 0.6. 

The transparency of the junction is chosen to be F = 
0.6A and the current is evaluated for U = 2.4A and 
U = 2.7 A. One observes that at relatively low bias volt- 
ages eV < 0.8A for U = 2.4A and eV < (0.9-=-0.95)A for 
U = 2.7 A the sub-gap current is essentially suppressed. 
For higher voltages the sub-gap current increases rather 
sharply, as a result of an interplay between Coulomb 
blockade and multiple Andreev reflections. The latter 
mechanism manifests itself in the occurrence of the sub- 
harmonic peaks in the differential conductance. Due to 
interaction, the positions of these peaks are shifted rel- 
ative to those in the noninteracting case eV = 2A/n, 
where n is the number of Andreev reflections. As can be 
seen in figure |^, increasing U results in a larger shift of 
peak positions. 

The behavior observed in figure ^ has a transparent 
physical interpretation which can be summarized as fol- 
lows. It is well known that in the absence of interac- 
tions the sub-gap conductance of superconducting junc- 
tions is determined by the process of multiple Andreev 
reflections (MAR) with the effective number of such re- 
flections n w 2A/eV. For the process with given n the 
charge {n+l)e is transferred between the electrodes. The 
relevant number n is obviously large at small voltages 
eV <^ A. Let us now turn on the interaction which 
strength is again characterized by some effective charg- 
ing energy Ec = e^/2C. Clearly, in this case MAR 
cycles with large n will be suppressed due to Coulomb 
effects much stronger than, say, one electron processes 
simply because in the case n 3> 1 the charge much larger 
than e is being transferred. E.g. at T = the sin- 
gle electron processes are blocked for eV < Ec, while 
MAR will be fully suppressed already at higher voltages 
eV < {n + l)Ec. This implies that for relatively small 
voltages the sub-gap conductance can be effectively sup- 
pressed even if the Coulomb energy Ec is small as com- 
pared to A. This is exactly what we observe in figure 



D 



FIG. 5. Same as figure |^ but for superconductors pair po- 
tential with p-wave symmetry. The parameters are U — 2.4 
(dotted curve ), U = 2.7 (sohd curve), e — —1.5 and F = 0.6. 

We observe that the sub-gap current for p-wave super- 
conductors is considerably larger than for s-wave ones, 
roughly by Imlx/Imax ~ 8. On the other hand, the ef- 
fect of the Coulomb repulsion U is rather similar. For 
U — 2.7A the current is suppressed compared to its value 
aiU — 2.4A. Beside the distinction of magnitudes, there 
is a novel additional structure in the I — V curves for p- 
wave superconductors which is related to the presence of 
a surface bound state. Comparing the results presented 
in Figs. 4 and 5 we observe that in the latter case, the 
current peaks at a certain bias voltage. This implies a 
negative differential conductance, which is the hallmark 
of resonant tunneling (contributed by the bound state) . 

Our analysis of the junctions formed by p-wave super- 
conductors can be straightforwardly extended to the case 
of d-wave pairing. The I—V curves and the sub-harmonic 
gap structure in junctions with d-wave superconductors 
in the absence of Coulomb interaction was recently stud- 
ied (see e.g. Ref. Q and other Refs. therein). Near zero 
bias the I — V curves | ]37| exhibit a current peak (equiv- 
alently, negative differential conductance) related to the 
presence of the mid-gap surface states. Notice, that in 
such systems the symmetry restricts the current, so that 
the contribution from the bound mid-gap states may van- 



ish if, for instance, one assumes 7^ ^ to be independent of 
Vx ■ As we have already argued before (see also Ref . ) , 
it might be essential to take the dependence of tunneling 
matrix elements on Vx into account already for point con- 
tacts. One can also consider the impurity model different 
from a point-like defect. Such a situation can be realized 
e.g. by artificially- induced defects Q . The spectroscopy 
of Bi2Sr2CaCu20s surfaces indicates that such defects 
appear to be more extended in STM imaging. In this 
case one can expect non-zero contribution from mid-gap 
level also in d-waves superconductors. Here, again, the 
electron-electron repulsion shifts the peak positions from 
their "noninteracting" values eV = 2A/n to higher volt- 
ages. It is quite likely that this interaction-induced shift 
was observed in the experiment [p9[. 



V. CONCLUSIONS 

In this paper the tunneling between two supercon- 
ductors or between a superconductor and normal metal 
through an Anderson-type quantum dot is investigated. 
Special attention is devoted to analyze the implications 
of the Coulomb repulsion between electrons in the dot on 
the tunneling process. The Andreev conductance for an 
SAN junction and the sub-gap current in an SAS junc- 
tion are calculated and elaborated upon. The theoretical 
treatment requires a combination of the Keldysh non- 
equilibrium Green function and path integral formalism 
and the dynamical mean field approximation. We derive 
general expressions for the effective action and the trans- 
port current through the system. These expressions are 
then employed in order to obtain a workable formula for 
the current. The latter is then calculated analytically 
and numerically for a certain set of energy parameters. 

The main results of the present research can be summa- 
rized as follows: 1) When one of the electrodes is a normal 
metal (an SAN junction) the gap symmetry structure is 
exhibited in the Andreev conductance. For p-wave su- 
perconductors, it shows a remarkable peak for voltages in 
the sub-gap region. For s-wave superconductors, on the 
other hand, the position of the peak is shifted towards 
the gap edge. It is further demonstrated that if the Hub- 
bard repulsive interaction increases the current peak at 
the gap edge is strongly suppressed. 2) The dynamics of 
tunneling between two superconductors (an SAS junc- 
tion) is more complicated. For s-wave superconductors 
the usual peaks in the conductance that originate from 
multiple Andreev reflections Q are shifted by interaction 
to higher values of V. They also suffer sizable suppres- 
sion as the Hubbard interaction strength increases. The 
sub-gap current in this case may describe the low energy 
channels in break junctions For p-wave supercon- 
ductors, the sub-gap current is much larger than in the 
s-wave case and the I — V characteristics exhibits a novel 
feature: the occurrence of mid-gap bound state results in 



a peak in the current, that is, a negative differential con- 
ductance. 
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VI. APPENDIX 



Below we will derive the result (62) within the frame- 
work of the formalism developed in the present paper. 
Consider a quantum dot between two s-wave supercon- 
ductors and assume that the interaction is negligibly 
small C/ — > 0. For the sake of simplicity we will also 
set Fi = Ffl = r. The result (^ can be expressed as a 
sum of two terms / = Iar + Iqp i where 



Iar — 



deTv[{NR)^'{gtr 



/-^-^ 

Jo 

12,-~^^21 / AT_\21/~R\12 I / M \21f~R\12-i 



-{NLY\g^Y' - {NrY'&IY' + {NLY'{g^^Y% (65) 
" Ah 



deTT[iNRY'C9t-9L 



\R\11 



{NLY\9R~~9Rr) 



(66) 



Here Iar is the sub-gap (Andreev reflection) contribu- 
tion to the averaged current while Iqp is defined by the 
excitations above the gap. In (A.1)-(A.2) we defined the 
Green-Kcldysh matrices g ~ ir^g with 



5«;f(e)^F«-(e)(e + r+A + r_A*), 



7R.A 



(67) 



VA2 - £2 



- A2 ' 

~9r,l{^) = [gld^) - ~91l{^)) tanh(e/2T). (69) 
We also defined 

Nl^ = {M^~gl^M^r'^, M^'^^ ^ {M"'^r\ (70) 

where the superscripts stand for the spin indices in 
Nambu space and Tr denotes the remaining trace over 
(discrete) energies which are scaled to A throughout this 
Appendix. 

Consider the limit of small voltages eV <C A. In this 
limit the sub-gap current Iar can be rewritten in the 
form 



Iai 



[9''{E^){{m\{M^r)\n) 



{n + l\{M^Y'\m)F^{E+) 
-{m\{M^Y^\n){n - l\{M'^Y"\m)F'^{E-)) 
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~g^{E-){{m\{MAY^)\n) 
{n\{M^f'\m)F^{E^) - {m\{M^y'\n) 
in\{M''f'\m)F''(Er,)) 
g''{E+){{m\{M^r\n) 
{Ti\{M^y^\m)F^{En) - {m\{M^f^\n) 
{n\iM^r\rr>,)F^iE^))]. (71) 

Here we denote E^ — eV{2n ± 1) and En — 2eVn. 
We also included r/2 into the definition of M and omit- 
ted terms non-diagonal in the spin indices because these 
terms are small in the limit eV <C A. At T ^ the 
summation over m is reduced to just one term with 
the maximum number mo determined by the condition: 

It is straightforward to evaluate the matrices 
{mQ\M'^'^\n) for sufhciently large F > A and eo — > 0. 
In this case M*'^ satisfy the following approximate equa- 
tions 

{m\{M"f'\mo)iElJ4-2) = 

+ (m - l|(M«)"|mo) + (m -f l|(A?«)"|mo), (72) 



- R 



{Em,F '{Em)E~ F iE° y^, 



(73) 



X R 



(m|(M )-|mo)«/4-2) = 



+ {m ~ l\{Mny^\mo) + (m + 1\{M R)^^\m„) (74) 

Similar equations can easily be derived for the two re- 
maining blocks. In the leading order in mo (this approx- 
imation is justified at small voltages 1^ ^ 0) at sub-gap 
energies (i?„ < 1, F^- — F^ = F) we obtain 



(n|(Al«)"|mo) = 



(-l)"(n + 1) 

(mo + 2)F(i;;^J ' 

(-l)«(n + !)£;„„ 



(n|(M^)2Vo) = - 



(mo + 2)^„F(i;„)' 
(-l)"(n + l)i;™o 



(mo + 2)^„i^(i;„)- 



(75) 
(76) 
(77) 



Substituting these matrix elements into eq. (A. 7) and 
performing a simple summation over n we arrive at the 
result (H). 
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